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Abstract 

The efficient simulation of models defined in terms of stochastic differential equa- 
tions (SDEs) depends critically on an efficient integration scheme. In this article, 
we investigate under which conditions the integration schemes for general SDEs can 
be derived using the Trotter expansion. It follows that, in the stochastic case, some 
care is required in splitting the stochastic generator. We test the Trotter integrators 
on an energy-conserving Brownian model and derive a new numerical scheme for 
dissipative particle dynamics. We find that the stochastic Trotter scheme provides 
a mathematically correct and easy-to-use method which should find wide applica- 
bility. 
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1 Introduction 



The study of mesoscopic particle models such as Brownian dynamics (BD)[1], 
Dissipative Particle Dynamics (DPD) [2,3], Smoothed Dissipative Particle Dy- 
namics (SDPD) [4] and the Voronoi fluid particle model [5,6] requires efficient 
integration methods that solve the appropriate stochastic equations of motion. 
In the past few years, several authors have considered improvements to the 
basic stochastic Euler schemes normally applied to these systems of equations, 
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particularly in the context of "conventional" DPD. Groot & Warren [7], Pag- 
onabarraga et al. [8] and Besold et al. [9] have reported various performance 
improvements to the basic schemes through the use of more sophisticated de- 
terministic solvers, for example those that have been successfully employed for 
deterministic dynamical systems [10] including molecular dynamics (MD) sim- 
ulations [11], such as the velocity and leapfrog Verlet algorithms. These tradi- 
tional deterministic integrators provide significant improvements on the basic 
Euler solver albeit, being deterministic schemes, their behaviour is completely 
uncontrolled from a theoretical point of view and their order of convergence 
is not clear. In fact, these solvers arbitrarily leave out terms which should 
appear in a correct stochastic expansion. More recently, alternative schemes 
have been devised resulting from proper stochastic expansions [12,13], and 
even from a Monte Carlo-based approach [14,15] where the fluctuations are 
introduced via a thermostat (the deterministic dynamics is still dependent on 
the integrator). 

A general method for deriving deterministic integrators is based on the Trot- 
ter expansion [1,16]. For Hamiltonian systems, these schemes preserve the 
symplectic structure of the dynamics and conserve the dynamical invariants, 
ensuring that the long time behaviour is correctly captured. In fact, if a dynam- 
ical invariant / exists then the discrete dynamics conserves exactly a virtual 
invariant /* which is bound to / up to second order in At [10]. An important 
feature of mesoscopic models is that they often recover a symplectic dynamics 
in some limit, an example being the DPD model for vanishing friction coeffi- 
cient. It may be important to account for this quasi- symplectic property of the 
SDEs in the integration scheme by assuring that in the same limit the scheme 
is symplectic as well [17]. 

Recently, a first order stochastic generalisation of the Trotter expansion has 
been rigorously proved [18,19]. In fact, for specific stochastic equations there 
exist schemes up to weak fourth order [20] or schemes corrected to reproduce 
more accurately the equilibrium distribution function [17]. The situation is less 
clear for a general SDE (such as Eq. (2) in Section 2), for which the application 
of the Trotter formula was overlooked in the literature, thereby generating 
some confusion in terms of how the Trotter formula can be used to split the 
stochastic equations. It is therefore useful to investigate the applicability of 
the Trotter formula in the most general case. This is of direct relevance for 
mesoscopic models which usually involve very large systems of SDEs. 

The Trotter formula has been applied to devise efficient integrators for sev- 
eral specific mesoscopic models but often its use is limited to splitting the 
propagator into several terms which are then integrated using standard nu- 
merical schemes. This approach would correctly produce the order of accuracy 
expected for the dynamics but potentially would affect adversely the conser- 
vation of the dynamical invariants or even detailed balance. Examples include 
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a numerical scheme suggested by a straightforward application of the Trot- 
ter rule to the Voronoi fluid particle model equations [21] which leads to time 
steps that are two orders of magnitude larger than the standard Euler scheme. 
In the context of the conventional DPD model, Shardlow [12,13] presented a 
new scheme, which splits the stochastic and deterministic parts following the 
Trotter rule, and then integrates the fluctuation-dissipation generators using 
the Bruenger et al. scheme [22] tailored onto the DPD equations. For Brownian 
dynamics, Ricci & Ciccotti [23] derived a numerical integrator based on the 
Trotter expansion which integrates the propagators by using the Suzuki for- 
mula [24] to transform the time-ordered exponential solution of the Brownian 
dynamics equations into more tractable simple exponentials. 



2 Stochastic Trotter schemes 

Let us consider first a deterministic dynamical system x(t) = C[x\. The formal 
solution of this system is x(t) = J2%Lo (^) p [ x ]( x o)(= e £ '[x](x )) as can be 
shown from the Taylor expansion around the initial condition xo. In general, 
the operator can be decomposed into simpler operators of the form C = J^f 1 £>%■ 
The Trotter formula (Strang [25]) provides a straightforward approximation 
to the time propagator 



where t = AtP, P is the number of time steps each of size At, and the 
ordering of the i,j indices is important. In the case that two operators A,B 
commute, i.e. [A, B] = AB — BA = 0, then the approximate Trotter formula 
is indeed exact because the equations e A+B = e A e B = e B e A are valid. Because 
the Trotter formula decomposes the dynamics over the time interval t into 
P steps, it provides a discrete algorithm for the solution of the dynamics of 
the system. Well known examples of the deterministic Trotter expansion are 
velocity and position Verlet schemes for molecular dynamics simulations [1]. 

In the stochastic case, we define a d dimensional stochastic process x t = 
{x\,...,xf) with associated stochastic differential equation (SDE) in the Ito 
interpretation 



where a fc (x t ) is the drift vector, 6 fej (x t ) is the diffusion matrix (d variables, m 
Wieners) and dW( the vector of independent increments of the j-th Wiener 
process. The mathematically equivalent Fokker-Planck equation (FPE) of Eq. 




(1) 



rn 



dx k t = a k {* t )dt + ]T b kj {* t )dWl 



(2) 
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Fig. 1. Diagram showing the derivation of SDE integrators by using the Trotter 
formula. Step (a) is the standard transformation from SDE to FPE formalism [26]. 
Step (b) decouples the FPE using the deterministic Trotter formula. Step (c) applies 
the reverse transformation from each decoupled FPE to the corresponding SDE in 
the order given in (b). 

(2) for the probability density p(x, t) is 

dtp = Hp\ (3) 

where T[p\ = - Ek £z (a h p) + \ £ M ^ (d M p) and d kl = £ j is the 
diffusion matrix. 

Following the diagram depicted in Fig. (1), we translate the starting stochastic 
equation (2) into the corresponding Fokker-Planck equation (3) which has 
formal solution p(x, t) = e :Ft [p](po)- The deterministic Trotter formula (1) can 
be applied to this formal solution by generally splitting the operator T = 
Furthermore, if Ti is a Fokker-Planck operator itself, this picture of 
evolving the probability density using the Trotter formula has a counterpart 
at the level of the SDE which would allow us to devise a numerical integrator. 
However, not all decompositions Ti have Fokker-Planck form and therefore an 
associated SDE. We then proceed by progressively splitting the terms in the 
starting SDE, i.e the drift vector a k and the matrix b k \ to verify Fokker-Planck 
form. 

The drift terms do not present any special problem: that is any splitting of 
the vector 

« fc = E<4 (4) 

a 

produces Fokker-Plank drift-like terms which can be easily integrated as with 
any standard ordinary differential equation (ODE). The diffusion operator 
demands more care. The matrix can be split into columns such as to give 
several systems of single noise equations, = b k i5 a j which are different 
from zero only in the column corresponding to noise a = j. By substituting 
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b kj = E Q b kj into the diffusive matrix d kl = J2 a ,/3 Ej b^b 1 ^ we obtain 

rf H = EE«' (5) 

a i 

which is split into several diffusive operators, because h a ■ bjg = 0,Vo; 7^ 
(3, i.e. the correlations between different diffusive dynamics are zero. In this 
procedure, we decouple the diffusive dynamics in terms of the subdynamics 
corresponding to each independent Wiener process. 

We are still left to integrate m single noise SDEs. We can try to decompose 
further each system of single noise SDEs into separate scalar SDEs. For each 
noise j, we set b k J = b kj 5 atk such that substituting in d kl we have 

d M = EEW. ( 6 ) 

a,P j 

which cannot be reduced to Fokker-Planck form for all terms. This means that 
we cannot split variables over terms of the same noise to derive the integrator. 
In fact, in order to apply the diagram of Fig. (I) and in particular step (c), we 
need to have all the terms in Fokker-Planck form to derive the corresponding 
SDEs. In principle, one could also try to separate the diffusion matrix d kl 
itself into several simpler matrices d kl = J2 a d& provided that each matrix 
is positive definite, but then the non-unique square-roots of the matrices 
have to be computed in order to recover the SDEs. Practically, this is very 
difficult in general. 

Finally, we must be able to compute the solution of the SDE corresponding to 
the i term Ti in order to write down the integration scheme. This is possible 
for simple SDEs, otherwise we can take advantage of the splitting between the 
drift and diffusion generators. The analytical solution of SDEs with zero drift 
is conveniently calculated in the Stratonovich interpretation for the stochastic 
integral (for a reference on Stratonovich integrals see [26]). In fact, the stan- 
dard rules of ordinary calculus apply and the SDEs are effectively integrated 
like ordinary differential equations by formally considering dW as dt. An Ito 
SDE like Eq. (2) is transformed into the equivalent Stratonovich form with 
the usual rules for the drift 

1 TO 

a k = a k -~Y,L j b kj (7) 

2 j=i 

where LP = Y%=i b h, ° g§7r an d the noise term is interpreted accordingly as 
dx k = a k (x t )dt + E™ 1 b kj (x t ) o dWl (see [26]). 

As the Trotter formula approximates the dynamics (3) of the probability dis- 
tribution p up to second order in time, we expect that at the SDE level the 
accuracy of the method is weak second-order [26], i.e. moments are accurate 
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to second order. Effectively, the proposed decomposition at the FPE level al- 
lows us to reduce the time-ordered exponential solution of SDE (2) in terms 
of simple exponentials up to second order provided that the generators for the 
same noise are not split. 



3 An energy-conserving Brownian model 

The oldest model for a stochastic system is the Langevin equation for a Brow- 
nian particle. In the one dimensional case, the SDE governing the velocity of 
the particle is dv = —vdt + (2T) 1 ^ 2 dW where we have selected units in which 
the mass of the particle and friction coefficient are unity and T is the di- 
mensionless bath temperature. This equation predicts an exponential decay of 
the velocity and, consequently, of the kinetic energy of the Brownian particle 
which goes into the fluid surrounding the particle. For illustrative purposes, 
we can construct an energy- conserving model in which we include the energy 
e of the fluid system, a Lagrangian reference system and a conservative force. 
We use the dimensionless equations in Stratonovich form 

dr = vdt 

dv = F(r)dt - vdt + (2ae) 1/2 o dW t , 

de = v 2 dt - (2ae) 1,2 v o dW t , (8) 

where F = — 9 Xjp- is the conservative force and a is a dimensionless heat 
capacity of the fluid. The above SDEs have as a dynamical invariant the total 
energy E = E = V(r) + y + e. Generalisations of the SDEs (8) to higher 
dimensions and multiple particles are indeed fundamental building-blocks of 
several mesoscopic models. 

In practice, it is not necessary to move to a Fokker-Planck description to derive 
the integration scheme. The derivation in section (2) shows that we can simply 
apply the Trotter formula (1) over the generators of the SDEs (8) provided 
that we do not split the stochastic generator for the same noise. The SDEs (8) 
is written in the form dx t = C[x]dt + «S[x] o dW t , where x = (r, v , e) and the 
deterministic and stochastic generators are respectively £ = £\+ £ 2 + £3 + £4 
and S = S\ + £2, 

£1 =vd/dr; C 2 = Fd/dv; £3 = —vd/dv; £ 4 = v 2 d/de; 

S 1 = (2ae) 1/2 d/dv; S 2 = -(2ae) 1/2 vd/de. (9) 

The generators S\ and c>2 cannot be split and integrated independently using 
the Trotter formula because they refer to the same noise. However, the solution 
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for S can be directly computed by applying standard calculus on the system 
of two equations dx. t = <S[x] o d\V t ] the solution is given by 



e SAw At [x] . e - C cos (v^A ^At + arccos(^/e, _ /y , ( , 

v — > sign(v)\^2C sin (-^/aAVt^ + arccos(y / e/C)^) , 

where sign(x) = 1 if re > and sign(x) = —1 if x < 0. Both variables are 
updated starting from the same initial values and C = e + v 2 /2 is computed 
before the update. The deterministic generators are easily integrated 



e £lA '[x] : r -> r + vAt; e £2A *[x] : u -> u + FAt, 

e £3A '[x] :v -*t;exp(-At); e £4A '[x] : e -> e + w 2 Al (11) 

The solutions of these differential equations can be nested following any given 
order to obtain different integration schemes. A possible numerical scheme is 

e S AW At/2 e U f e C 3 M e C 2 f e d At e C 2 f e C 3 M g £ 4 M g S AW' At/2 ^ ^ 

where AW' At , 2 and AWaj/2 are two random numbers drawn from a zero mean 



normal distribution with standard deviation y At/2. We note that the stochas- 
tic propagator of this scheme conserves energy exactly (for any time step size), 
therefore the conservation of energy depends only on the approximation in- 
troduced in the deterministic part. 

As already stated, it is not possible to decompose the stochastic generator S 
into two independent stochastic scalar equations using the Trotter formula. 
Unfortunately, this approach is what would follow if one was to apply naively 
the Trotter formula to SDE (8). The resulting scheme would not be second 
order and would conserve energy poorly. For instance, this is the case for the 
scheme 

e Si AW At/2 e S 2 AW At/2 e U % e C s % e C 2 % e CiAt e C 2 % e <C 3 ^ e U % gSa AW^ t/2 ^A^^ 

(13) 

where the stochastic propagators are 
e SlAWAt [x\ :v^v + V2^~eAW At , 

e S2Alfil [x] : e -> (^i - Vtov/2AW M ) 2 . (14) 



Interesting, there is a possibility to apply a Trotter-like rule to devise second 
order weak integrators even for the decomposition S = Si +S2- To do this the 
noises have to be advanced by Av ^ At = (weak) AW At/ 4, where by = (weak) 
we mean that moments of both sides are equal to second order. Note that for 
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Fig. 2. Average of the maximum relative error of the energy for the SDE (8) over 
10 independent runs up to t = 1 for the Trotter schemes (12) circles, (15) squares 
and the incorrect Trotter scheme (13) triangles. The deterministic Trotter scheme 
for Eq. (8) with a = 0, vq = 1, eo = 1/2 is plotted with dotted lines for reference. 

the Trotter expansion it should be AW At/2 = W t+ At/2 — W t . The scheme is 
written as 

(15) 

where we use the same realization of the noise AW^t/A- The second order 
weak convergence can be verified by a direct comparison with a second or- 
der stochastic expansion and intuitively understood by formally considering 
AW as At. We stress that the resulting scheme does not correspond to a 
stochastic Trotter expansion, but rather to a second order approximation of 
the propagator. This method provides a way to write an integration scheme 
even in cases where it is impractical to compute the solution of the generator 
S altogether. However, wherever possible, this approach should be avoided 
or limited to the smallest generator because the resulting integration scheme 
may loose important structural features of the dynamics (as in the example 
of SDEs (8)). 

We validated numerically the integration schemes (12) and (15) as well as 
the incorrect one (13). The simulations were run using the bistable potential 
V(r) = f3{r A — 2r 2 ) with a — 1, (5 — 1 and initial conditions r = 0, v — 
and eo = 1. The average relative error for the total energy AE/E for dif- 
ferent time step lengths At is shown in Fig. 2. The error is computed by 
averaging the maximum error reached by t = 1 over 10 independent runs. The 
stochastic- Trotter scheme (12) conserves the energy with the same accuracy 
as the deterministic Trotter scheme (computed using a = 0). The scheme (15) 
is consistent with first order accuracy (it is second order for single time step 
error), while the incorrect scheme (13) does not conserve energy with first 
order accuracy. Note that the order for the cumulative error is one less than 
the single time step error. Clearly, the energy conservation performance of 
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the Trotter scheme (12) is a direct consequence of the exact integration of its 
stochastic component which is impossible to achieve by other general schemes. 



4 A Trotter integration scheme for dissipative particle dynamics 



We now apply the stochastic Trotter expansion to the equations of dissipative 
particle dynamics. The DPD model consists of a set of N particles moving in 
continuous space. Each particle k is defined by its position r k and its momen- 
tum p k and mass m. The dynamics is specified by a set of Langevin equations 
very similar to the molecular dynamics equations, but where in addition to 
the conservative forces there are dissipative and fluctuating forces as well 



dr k = p k /mdt, 

N 

dp k = J2 e M [akiF c {r k i)dt - ~f/muj D (r k i)(e k i ■ p k i)dt + au^r^dW^ ,(16) 

where F c (r) is the conservative pair interaction force weighted by positive and 
symmetric parameters a k i, r k i — r k ~ r i is the distance between the particle 
k and particle /, r k i its length and e k i = r k i/r k i. The weight functions ud,^r 
usually have finite range r c and are related by u>D(r k i) = u R (r k i) in order to 
satisfy detailed balance. This condition ensures that the equilibrium state is 

2 

Gibbsian and sets the value of its temperature to T = ^k^- ^ typical selection 
is LU R (r k i) = io(r k i) with 



1 — — r < r, 



c 



r > r c . 



(17) 



The conservative force F c (r k i) = — 9V qI^ is usually chosen to be of the form 
F c (r k i) = w(r k i). 

The generator of DPD equations (16) is £ = + £ fc>J#fc (v kl + S kl ) , 

where 



C k r = p k /md/dr k ; S kl = au R (r kl )e kl d/dp k ; 

V kl = a M F c (r k i)e kl d/dp k - >y/muj D (r k i)(e k i ■ p kl )e k id/dp k . (18) 

In the DPD model the momentum is conserved because the forces between in- 
teracting particles k and I satisfy Newton's third law. We split the DPD equa- 
tions in order to satisfy this requirement. The conservative and fluctuation- 
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dissipation generators for the pair interaction k, I give 

rfx = (V kl + V lk ) [x]dt + (S kl + S lk ) [x]dW& (19) 

where x = (ri, r N , pi, pjv). The solution is computed by noting that 
dpk + dpi = and dpk = \dpui where dpki = dpk — dpi. The equation for dpki 
can be solved for the component of the radial direction because from the form 
of the SDEs (16) it follows that dpki = d(p k i -e^e^. Let us call p\ x = pui -e^; 
then we have an Ornstein-Uhlenbeck process 

dp e kl = Adt - Bp e kl dt + CdWl l: (20) 

where A = 2akiF c (r k i), B = 2^/mujD and C = 2<jur, which has analytical 
solution [26] 

Pti(t) = e- BAt pti(t ) +Af e B ^ds + C f e B ^dW s , (21) 

Jt J to 

where At — t — to, to being the initial time. The solution (21) of the Ornstein- 
Uhlenbeck process requires the generation of coloured noise based on a nu- 
merical scheme itself [27]. In fact, the stochastic process p%i(t) has stationary 
correlation function for t, s — > oo with finite \t — s\ given by 

< PtiWM >= ^ + ^ exp(-5|f - S |). (22) 

A version of the method to generate coloured noise [27] adapted to Eq. (21) 
results in the scheme 



Apll = ( m . etl _ *£k) _ t ) + *^f" (23) 



where (, kl = £ lk are normal distributed with zero mean and variance one 
(7V(0,1)) and Ap& =pfc(f) -pfcft,) . 

The propagator K kl for p k and p ; is then given by 

£&[x] : (p fc , p,) -> (p fc + ^Ap|,e«, Pi - ^Ap|,ejw) . (24) 

The remaining position update is given by 

e^ At [ x ] : r fc -> r fc + Pfe /mAt. (25) 
We note that commutes with therefore we can use the exact formula 

e ^^ At = n^=ie £ * At . 
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The DPD scheme is finally given by the following Trotter integrator 

N N 1 

x(t+At)^ n n^ At n ^u*^- ^ 

k=l,l=l k=l k=N,l=N 

In practice the integration algorithm consists of the following steps: for the 
interaction pairs k,l update the momentum half timestep according to the 
propagator (24), where £ kl = t; lk are drawn from a normal distribution with 
zero mean and variance one; iterate over particles k updating the position 
according to (25); finally, update pairs k,l in reverse order again using the 
propagator (24) but with new noise £' kl . This algorithm requires the calculation 
of the pair-list only once per iteration and has the same complexity as a simple 
DPD velocity- Verlet scheme (DPD-VV [7]). 

We test this integration scheme using the open-source code mydpd [28] written 
in simple C++ and implementing the DPD models described here with peri- 
odic boundary conditions. The simulations are run with N = 4000 particles, 
aw = 25, 7 = 4.5, a — 3, m — 1, r c — 1 in a three dimensional periodic box 
(L, L, L) with L = 10. These settings give a particle density p = 4 and equilib- 
rium temperature /c#T = 1. In our implementation, the computational cost of 
each scheme averaged over several iterations indicates that the Trotter scheme 
is 60% more costly than the simple DPD-VV but 10% faster than the Shardlow 
SI scheme (which costs almost twice than DPD-VV). The equilibrium temper- 
ature for the DPD-Trotter scheme of Eq. (26), DPD-VV [7] and Shardlow [12] 
schemes is reported in Table 1. The DPD-Trotter scheme recovers the equi- 
librium temperature better than DPD-VV, and as accurately as Shardlow's 
scheme. This difference depends on the implicit scheme used by Shardlow for 
the integration of the pair interaction. In our case, we have used an exact 
integration Eq. (21) which, however, requires the generation of coloured noise 
[27] which is by itself a numerical scheme. Considering the accuracy of the 
equilibrium temperature and the computational cost, both DPD-Trotter and 
Shardlow schemes are integrators of comparable performance for the DPD 
equations. A more detailed study of the equilibrium properties of the fluid is 
necessary to assess the accuracy in reproducing the equilibrium distribution 
and other statistical properties. 



5 Conclusions 

The stochastic Trotter schemes can provide efficient integrators for stochastic 
models with dynamical invariants by fully taking into account the underlying 
stochastic character. The stochastic Trotter formula can be applied to any 
model based on SDEs and should find wide applicability provided that some 
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Table 1 

Equilibrium temperature for the DPD- Trotter, Shardlow and DPD-VV schemes for 
different time steps. The average of the kinetic temperature < ksT > is computed 
over a simulation of duration t = 1000. The standard deviation of the estimates, 
computed by block-averaging, is less than ±5 x 10~ 4 . 



At DPD-Trc 


>tter (scheme Eq. (26)) ! 


Shardlow [12 


] DPD-VV [7] 


0.05 


1.0136 


1.0138 


1.0411 


0.02 


1.0020 


1.0018 


1.0097 


0.01 


1.0007 


1.0005 


1.0043 



care is used to decouple the stochastic dynamics for the same noise. These 
types of stochastic schemes offer the flexibility to easily tailor the integrator 
to the specific model, thereby integrating exactly important parts of the dy- 
namics. This stochastic Trotter scheme is a second order weak scheme, but, 
more important, in our examples it provides very good conservation of the 
dynamical invariants. 
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